MULTIPLE QUENCHING SOLUTIONS OF A FOURTH ORDER 
PARABOLIC PDE WITH A SINGULAR NONLINEARITY 
MODELLING A MEMS CAPACITOR. 



A. E. LINDSAY * AND J. LEGAt 

Abstract. Finite time singularity formation in a fourth order nonlinear parabolic partial dif- 
ferential equation (PDE) is analyzed. The PDE is a variant of a ubiquitous model found in the 
field of Micro-Electro Mechanical Systems (MEMS) and is studied on a one-dimensional (ID) strip 
and the unit disc. The solution itself remains continuous at the point of singularity while its higher 
derivatives diverge, a phenomenon known as quenching. For certain parameter regimes it is shown 
numerically that the singularity will form at multiple isolated points in the ID strip case and along 
a ring of points in the radially symmetric 2D case. The location of these touchdown points is ac- 
curately predicted by means of asymptotic expansions. The solution itself is shown to converge to 
a stable self-similar profile at the singularity point. Analytical calculations are verified by use of 
adaptive numerical methods which take advantage of symmetries exhibited by the underlying PDE 
to accurately resolve solutions very close to the singularity. 
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1. Introduction. Micro-Electromechanical Systems (MEMS) combine electron- 
ics with micro-size mechanical devices to design various types of microscopic machin- 
ery (cf. [29] ) - A key component of many MEMS is the simple capacitor shown in 
Fig. 11.11 The upper part of this device consists of a thin deformable elastic plate 
that is held clamped along its boundary, and which lies above a fixed ground plate. 
When a voltage V is applied between the plates, the upper surface can exhibit a 
significant deflection towards the lower ground plate. When the applied voltage V 
exceeds a critical value V*, known as the pull-in voltage, the deflecting surface can 
make contact with the ground plate. This phenomenon, known as touchdown, will 
compromise the usefulness of some devices but is essential for the operation of others 
(e.g. switches and valves). Capturing and quantifying this phenomenon is a topic of 
some mathematical interest and is the subject of this paper. 



Free or supported boundary 



Elastic plate at potential V 




Fixed ground plate 



Fig. 1.1. Schematic plot of the MEMS capacitor (reproduced from \27j ) with a 
deformable elastic upper surface that deflects towards the fixed lower surface under an 
applied voltage. 
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A canonical model, originally proposed in |29j . suggests the dimcnsionless deflection 
u(x, i) of a device occupying a bounded region O C R 2 satisfies the fourth-order 
problem 

\f(x) u = 0, d n u = x^dVt. 

ut = — A it + 6Au — * n > , (1.1) 
(l + u) 2 « = 0, t = xen 

Here, the positive constant 5 represents the relative effects of tension and rigidity 
on the deflecting plate, and A > represents the ratio of electric forces to elastic 
forces in the system, and is directly proportional to the square of the voltage V 
applied to the upper plate. The function f{x) G C a (Vt) for a G (0,1), represents 
possible heterogeneities in the deflecting surface's dielectric profile while the boundary 
conditions in (jl.lj) assume that the upper plate is in a clamped state along its rim. 
The model (jl.ip was derived in |29| from a narrow-gap asymptotic analysis. 

The second order equivalent of (jl.ip 

Xf(x) u = 0, xedfl. 

(i + u) 2 u = o, t = o xen 

has been the subject of extensive study recently and there are now many established 
results regarding the behaviour of solutions, both dynamic and steady (c.f. [7] and 
the references therein for a thorough account). In particular it is known that there 
exists a A* > such that whenever A > A* and inf^ / > 0, the device touches down 
in finite time, i.e. ||1 + u(-,t)||i n f — > + as t — > t~ . Lower and upper bounds have 
been established on the touchdown time t c of (|1.2j) and it is known that if touchdown 
occurs at an isolated x c G il, then f(x c ) ^ 0. Additionally, a refined asymptotic 
study of the touchdown profile was performed in |15j where it was shown that the 
quenching solution is not exactly self-similar and has asymptotic form 

U -> -1 + [3/(x c )A(t c - t)] 1/3 ( 1 - — j- rr + (X ~ rr + ■ ■ ■) 

1 1 >y n \ 2|log(* c — 1)| 4(t c - i)|log(t c - i)l / 

(1.3) 

where x c G f2 and t c > are the touchdown location and time respectively. In 
addition, when f(x) is a constant and O = [—1,1], the unique touchdown point is 
Xc. = 0. 



In contrast to the second order problem (|1.2|) . very much less is known about the 
fourth order problem (|l.ip . partly due to the lack of a maximum principle. In the 
absence of the tension term (6 = 0) and with f(x) = 1, equilibrium solutions of (|1. 11) 
were studied in [T^ and the existence of a pull in voltage A* was demonstrated for ft 
a radially symmetric ball. The maximal branch of equilibrium solutions to (jl.ll) , i. e. 
those solutions with largest L 2 norm for any sufficiently small A, were constructed 
in the limit as u —> — 1 + in one and two dimensions in [3SJ [55] ■ Under the relaxed 
Navier boundary conditions u = Au = on 951, a maximum principle is available 
and theoretical results regarding the existence and uniqueness of solutions are more 
tractable IT71I2T1. 



Literature on the dynamics of fourth order MEMS equations is particularly sparse. 
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The work of |14j concerning the wave equation 

HW tt +w t -Aw + BA 2 w = in SI x (0, T] 

w = Aw = on dttx(0,T) (1.4) 

w(x, 0) = wq(x), Wt (x, 0) = wi (x) in f2 

appears to be the first contribution to the topic in which it is shown that (|1.4p touches 
down in finite time for A > A*. 

In the present work radially symmetric dynamical solutions of the fourth order MEMS 
problem 

u t = -e 2 A 2 u - - — ~ — ~T , iGfi, u(x, 0) = 0, x 6 f2 (1.5a) 
(1 + u)^ 

are considered for domains 

(Strip) : fi = [-1, 1]; (Unit Disc) : O = {.t 2 + y 2 < 1} (1.5b) 
and boundary conditions 

(Clamped) : u = 0, 9„u = x £ <9f2; (Navier) : u = 0, Am = x 6 90. 

(1.5c) 

The particular form of this equation is obtained from (jl.l[) by setting f(x) = 1, 
neglecting the tension term Am (5 = 0), taking At as a new time variable, and defining 
A = e~ 2 . The consideration of radially symmetric solutions of (|1.5aj) on the strip and 
unit disc geometries effectively focusses attention on the PDE 



Ut = -e 2 



„„ 2(jv-i) ,„ i^-i „ i^-i ; 

u H u = — u H 5 — u 



(1 + u) 2 



(1.6) 



for iV = 1 ( Strip ) and iV = 2 (Unit Disc). 

The paper begins with some proofs confirming that ()1.5|) exhibits the pull-in insta- 
bility, i.e. there is a number e* > such that when e < e* , (jl.5p has no equilibrium 
solutions and will touchdown to u = — 1 in finite time. In a moving mesh PDE 
method (MMPDE) is employed together with an adaptive time stepping scheme to 
accurately resolve the solution of (|1.5j) very close to touchdown. While touchdown 
occurs at the origin for certain parameter regimes as in the second order equivalent, it 
is observed that for e below some threshold £ c , equation (| 1 . 5|) may touchdown at two 
separate isolated points in the strip case and, under radially symmetric constraints, 
along a ring of points in the unit disc case. Moreover, it is observed that the location 
of the touchdown set has a dependence on e that can be analyzed. While multiple 
touchdown has been observed previously when tailored dielectric profiles f{x) were 
considered, here the device is uniform ( f(x) = 1 ) and the location of touchdown 
can be parameterized through e = A -1 / 2 . This may potentially allow MEMS devices 
to perform more exotic tasks or simply extend their lives by spreading wear over a 
larger area. 

In 21 the location of touchdown for (|1.5|) is analyzed by means of asymptotic expan- 
sions which predict that that for the strip case, the two touchdown points are 



1 - £ 1/2 /(i c ) 1/4 [% + f(t c )m + f 2 (t c )V2 + ■■■]], f(t) = 1 - (1 - 3^ 3 

(1.7) 
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while for the unit disc radially symmetric touchdown occurs on a ring with radius 



1 



(1.8) 



where r/o, r]i , ?72 and rji , rji are numerically determined constants whose values depend 
on the boundary conditions applied (|1.5b|) . Note that these asymptotic predictions 
are in terms of the touchdown time t c and are valid for e < e c . In order to estimate the 
values of x c and r c , a numerical approximation of t c is required. These formulae are 
shown to agree well with full numerics, particularly when e«1. The limiting profile 
of (|1.5[) as inf^go u(x, t) — > —1 is also constructed. In contrast to the quenching 
profile (|1.3|) of the second order problem (|1.2j) . it is observed that (|1.5j) exhibits a 
self-similar quenching profile which finalizes to 



u(x, t) 



1 



co 



—Trr 



4/3 



MS 



(1.9) 



where the parameter cq is determined numerically and has value cq = 0.9060 for both 
the strip case and touchdown away from the origin in the radially symmetric unit disc 
case. In the unit disc geometry with touchdown at the origin, the numerically obtained 
value is cq — 0.7265. The stability of this profile is determined and convergence of 
the numerical solution of (|1.5|) to the self-similar profile (|1.9[) is verified in each case. 

2. Preliminary Results. In this section two preliminary results are established. 
The first result demonstrates that for e small enough, (jl.5|) has no equilibrium solu- 
tion. The second result proves that when no equilibrium solutions exist for (|1.5[) . the 
solution will touchdown, i.e. reach u(x,t) = — 1, at some point in space in some finite 
time. These results rely on a positive eigenpair (<^0jW)) of the problem 



A 4> = [j,(f>, x g fi; 
for the strip and unit disc geometries and the boundary conditions 



(2.1a) 



(Clamped) : = 0, d n <f> = x G <90; (Navier) : = 0, A</> = x G 90 

(2.1b) 

In the case of clamped boundary conditions, it is well known that for general two- 
dimensional geometries, the principal eigenfunction of (|2.ip need not be of one sign. 
Two well known cases are that of the square [5] and annulus [B]. However, if only 
the strip and the unit disc are considered, then (|2.1[) docs admit a strictly one signed 
principal eigenfunction together with a positive eigenvalue. A brief calculation shows 
that the eigenfunctions for the clamped strip satisfy 



4> = C 



sin£(a; — 1) — sinh^(a; — 1) 



sin 2£ — sinh 2£ 



cos 2£ — cosh 2£ 



where £ = /i 1 / 4 and 



cos 2£ cosh 2^ = 1. 



[cos£(.t — 1) — cosh£(a; — 1)] 
(2.2a) 

(2.2b) 



For the clamped unit disc, the eigenfunctions are 



= C 



MOW = Wo(0- 



(2.2c) 
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where again £ = /i 1 / 4 . In equations (|2.2[) . the constant C is fixed by normalization. 
The case of Navier boundary conditions for this eigenvalue problem were considered 
in |14j where it was shown that (</>OjMo) = (<Pn, ^n) for 

A0 n + A n <fo = 0, x G fi; <fo = 0. i e dft (2.3) 

is a positive eigenpair of (|2.1a[) . The maximum principle guarantees the positivity of 
the principal eigenpair of (|2.3j) for any f2 C [15j . The principal eigenfunction for 
the strip and unit disc geometries under Navier boundary conditions are therefore. 

(Strip): Mo = ^ <£o = Can. (|(x - 1)) (2.4a) 
(Unit Disc) : fj, = z\ cj) = CJ (z r), (2.4b) 

where C is a normalization constant and in (|2.4b[) zq is the first root of Jq(zq) = 0. 

The following theorems show that for e small enough, equation (|1.5[) admits no equi- 
librium solutions and will touchdown to u = — 1 in finite time. The proof techniques 
involved have been employed previously in |19jl27j and rely on a positive eigenfunction 
of (|2.ip . Therefore, in the case of clamped boundary conditions for (| 1 . 5|) . the result 
is limited to the strip and unit disc geometries. 

Theorem 1 : (c.f. JZ21 \F7p There exists a real < e* < oo such that for < e < e* , 
equation (|1.5|) has no equilibrium solutions when considered on the strip or unit ball 
with clamped conditions and any Q, C R 2 for Navier conditions. In addition e* > e — 
a/27/4/xo where (</>o,Mo) * s a positive eigenpair of (|2.1[) . 

Proof : Take (4>o, /iq) to be an eigenpair of (|2.1[) with <f>Q > and > 0. Multiplying 
the equilibrium equation of (|1.5[) (i.e. u t = ) by 0o and integrating gives 

0o W + n+^a ) dx = ( 2 - 5 ) 

Clearly (|2.5|) cannot hold when the integrand is strictly positive which occurs when 
the inequality 

e 2 iM)U + 1 > (2.6) 
(1 + it) 

is satisfied on f2. This implies that e* is finite. The equality e 2 fiQU = —(1 + u)~ 2 
has exactly one solution when e 2 /io = 27/4, and no solutions when e 2 /io < 27/4. 
Therefore, whenever e 2 Mo < 27/4, (|2.6[) holds and (|1.5j) certainly has no equilibrium 
solutions. Moreover, the smallest positive e such that (|1.5[) has an equilibrium solu- 
tion, e* satisfies 



Numerical values of Mo, satisfying the first positive solutions of (|2.2[) . e and e* are 
given in Table. [2~T1 under both clamped and Navier boundary conditions. I 



The following theorem shows that for e < e, when an equilibrium solution to (|1.5[) is 
not present, touchdown occurs in finite time. 
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Theorem 2 : Suppose that e < e — y / 27/4/x , then the solution of (|1.5|) reaches 
u = — 1 in some finite time t c when considered on the strip or unit ball with clamped 
boundary conditions, or on any bounded Q C K 2 under Navier boundary conditions. 

Proof : The proof follows Theorem 3.1 of [19] and relies on the existence of a positive 
eigenfunction </>o of ()2.1[) . Let (j>o be normalized by the condition J Q (j)odx = 1. By 
multiplying (|1.5a[) by <f>o and integrating by parts, the equality 

d 



(j) udx = -e 2 fj, / 4> udx- / — — — dx (2.8) 
dt Jn Jn Jn ( l + u ) 2 

is obtained. Defining E(t) = J n ipgudx where 22(0) = implies that 

dE o „ f d>Q 9 „ 1 

-e 2 /ioJ5 - / r-, ^ dx ^ -e IM)E - 



dt ^' J n (l + u) 2 ~ ^ (l + E)^ 

where in the last step Jensen's inequality has been applied. Standard comparison 
principles show that E(t) < F(t) where F(t) satisfies 

dF 1 

¥ = -^-(iW' F(0) = ' (2 ' 9) 

Equation (|2.9p is separable and so it is solved to show the touchdown time i, for F(t) 
at which F(t) = —1 satisfies 

i= [ (s 2 fi s+ ., 1 yr ) ds. (2.10) 



-l 



(1 + 8)* 



The touchdown time for F(t) is finite when this integral converges which occurs when 
e < e = y / 27/Afi . Finally, since 

E(t) = \ (f>oudx > inf u (f>a dx = inf u. 
Jn xen Jn xen 

it follows that 

inf u < E(t) < F(t) (2.11) 
x^n 

so that if i from (|2.10[) is finite, then the touchdown time of (|1.5p . t c must also be 
finite. Therefore when e <= ^/27/4/xo, t c < i where i is given in p.lOp . In the limit 
as e — > + , equation (|2.10p has expansion 

t=l + £ -W + °(^ ( 2 - 12 ) 



This provides the asymptotic upper bound on the touchdown time t c of (jl.5 

3 + 30 

in the limit as £ —> + . 



t c<l + £ -^+0(£ 4 ) (2.13) 



The preceding analysis demonstrates the presence of the ubiquitous pull-in instability 
for (|1.5[) for general geometries when the boundary conditions are Navier and for the 
strip and unit disc domains when clamped boundary conditions are applied. It is 
an open and challenging problem to prove that (|1.5p exhibits the pull-in instability 
for general geometries il C when clamped boundary conditions are applied. A 
useful byproduct of the analysis presented here is the estimators on the critical pull-in 
voltage e* for each of the geometries and boundary conditions considered as collated 
in Tabic O 
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Mo 


£ 


e* 


Navier B.C.s 
Clamped B.C.s 


Strip Unit Disc 
6.0881 33.4452 
31.2852 104.3631 


Strip Unit Disc 
1.0530 0.4492 
0.4645 0.2543 


Strip Unit Disc 
1.0771 0.4695 
0.4778 0.2683 



Table 2.1 

Numerical values of the principal eigenvalue fio from 112.21 1. s from Theorem 1,2 
and e* under clamped and Navier boundary conditions. The values of e* were calculated 
numerically in \27j as a saddle-node bifurcation point of equilibrium solutions to 11.51 . 



3. Numerics. In order to obtain accurate numerical representations of (|1.5p 
close to touchdown, a method which can resolve the rapidly changing spatially lo- 
calized and temporal features of the equation seems warranted. To facilitate this, 
the r-adaptive moving mesh scheme MOVCOL4 of [31] together with the adaptive 
time stepping scheme of [5] are implemented. Both schemes take advantage of the 
underlying invariance of equation (jl.5l) to the transformation 

t->at, (1 + u) -> a 1/3 (l + u), x->a 1/4 x. (3.1) 

A brief overview of the method is now provided, for more details see [31 [3T|. The 
physical domain is approximated by the grid 

.t < xi(t) < ■ ■ ■ < xjsr(t) < xn+i, (3.2a) 

the node points of which are evolved with the equation 

-"/X m = (M(X)X s ) s . (3.2b) 

Here 7 is a small parameter which controls the relaxation timescale to the equidis- 
tribution profile, M(X) is known as the monitor function and Xj(t) = X(iA£,t) is 
a map between the physical domain and a computational domain fi c = [0, 1] with 
coordinate £ € [0, 1]. In calculations, the value 7 = 10~ 4 was used and the boundary 
conditions Xq = X^ + i = were applied. The monitor function 

was selected which provides a balance between grid points in the region where M is 
large (e.g. where ||1 + u||i n f is small) and also in regions where the solution is not 
changing rapidly but modest resolution is still required so that iterative procedures 
converge. Importantly, with this choice of M(u), equation p.2bp retains the sym- 
metry (|3.1[) of the underlying equation. Spatial discretization was effected by a 7 th 
order polynomial collocation procedure with evaluation at four Gauss points in each 
subinterval (c.f. Appendix |A"| . After accounting for boundary conditions, this results 
in a system of A(N + 1) equations for the solution and its first three derivatives at 
each node point. The mesh equation (I3.2b[) is discretized as follows 

AVi - 2Xi + X l+1 M l+ i (X i+1 - Xi) - M i _i(X t -Xi-i) 

- ^ ^ = — Ae (3 - 4a) 
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where 

M i+ ,_ = M(* +1 ) + M(*) (34b) 
2 2 

The integral term of (|3.3j) is evaluated by the trapezoid rule on the subintervals defined 
by the XiS. The efficient simulation of the PDE close to singularity necessitates the 
use of temporal adaptivity. The underlying symmetry of the problem (|3.1[) provides 
indication on how the time stepping should be adjusted according to the solution 
magnitude and motivates the introduction of a computational time coordinate 

~T~ = 9{u) g(u) = — ( \, r , (3.5) 

dr mi xe n\\M(u)\\ 

where again (|3.5p retains the underlying symmetry (|3.1[) of the underlying problem. 
The discretized main equation (|1.5[) and equations for the mesh ()3.2b[) arc written in 
terms of the computational time r and solved simultaneously as a DAE of form 

= M(y,T)y T -f(y,T), y = (t(r), u, X) T . (3.6) 

Here u G M 4 ( Ar + 2 ) is a vector containing the nodal values of the solution and its first 
three derivatives while X £ l w+2 is the vector of grid points. The square mass matrix 
Ai is of size 5(N + 2) + 2 and has entries filled with the discretizations of (|1.5[) and 
(|3.3j) while f € R 5 ( w+2 )+ 2 represents the discretized right hand sides. The resulting 
equations are solved in MATLAB with the routine ode23t. 

In Fig. 13. li the three solution regimes for (|1.5|) on the strip under clamped boundary 
conditions are observed. When e > e* the beam attains a steady equilibrium deflection 



and does not touchdown (c.f. Fig. 3.1(a) I. The second solution regime lies in the 



parameter range e c < e < e* whereby the solution touches down in finite time at the 



origin only, as displayed in Fig. 3.1(e) The simulation is halted when inf^gn ||1 + 
u(x, t)\\ reaches a specified proximity to u = —1. In the case N = 1 with e = 0.2, the 
solution can be followed to u(0) = -0.99999 with t c - 1 = O(10" 17 ). In the case of 
multiple touchdown points symmetric about the origin, the solution can be followed 
to infsefi = -0.999 where t c -t = O(10" 10 ). When multiple touchdown 

points are present it is more challenging to integrate (|1.5|) very close to touchdown 
as grid points will tend to coalesce on one of the two touchdown points thereby 
hindering convergence at the other. On the figures displaying numerical solutions, the 
grid points are indicated on the curve as crosses and are observed to coalesce on the 



singularity point as t — s- t c (c.f. Fig. 3.1(d) I. In the third parameter regime < e < e c , 



touchdown occurs in finite time at two isolated points symmetric about the origin (c.f 



Fig. 3.1(e) I with the location of touchdown as a function of e indicated in Fig. 3.1(b) 
The border of the one and two point touchdown regimes is approximately e c ~ 0.066. 

In the radially symmetric unit disc case, touchdown occurs at the origin when e c < 
e < e* and on an inner ring of points when e < e c ~ 0.075. 

A possible interpretation for this behavior is that u = —1 is an attractor of the system 
and that the location of touchdown is governed by the critical points of the deflection 
u{x,t) as the solution enters the basin of attraction for u = —1. This would suggest 
that the source of the multiple touchdown points lies in the dynamics of (jl.5p for 
small t. 



4. Asymptotics. 
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(a) e = 0.5 (b) Touchdown Points 




(c) e = 0.2 with N = 16 (d) e = 0.2 with N = 16 Zoomed. 




(e) e = 0.02 and N = 24. (f) e = 0.02 and N = 24 Zoomed. 



Fig. 3.1. The above figures relate to numerical solutions of HI, 51 for the strip 
domain with clamped boundary conditions. The mesh points are indicated on solutions 
with small crosses so that their dynamics can be observed. In panel (a), solutions are 
shown for e = 0.5 > e* so that touchdown does not occur and a steady state deflection 
is approached. Panel (b) displays the relationship between touchdown location(s) and 
the value of e. The critical value e = e c , below which touchdown occurs at two points, 
is approximately e c = 0.066. In panels (a) and (c)-(f), solutions are increasing in time 
from top to bottom. Panel (c) shows solutions for e = 0.2 < e* and touchdown is 
observed at the origin around time t = 0.3833. In Panel (d) a zoom in of the touchdown 
region is displayed which shows the refinement of the mesh in this area. In panel (e), 
solutions are shown for e = 0.02 where touchdown is observed at two separate points, 
symmetric about the origin around t = 0.3240. Panel (f) displays a zoom in of the 
positive touchdown region for Panel (e) and again the refinement of the mesh in this 
region is apparent. 
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4.1. Small time asymptotics. In this section an analysis of the biharmonic 
MEMS equation 

u t = -e 2 A 2 u- - 1 x9 , xeQ; u(x,0)=0, xefl (4.1) 

(i + uy 

is performed in the small time regime t — > + for strip and disc domains (|1.5b| and 
boundary conditions (|1.5cp . In this regime the deflection of the beam is small which 
allows the (1 + u)~ 2 term to be linearized and in this way its influence can be thought 
of, to leading order, as a uniform forcing term of unit strength. 

In a region away from the boundary where the e 2 A 2 u is negligible for e <C 1, the 
leading order solution satisfies u{x, t) = u(t) where 

ut = ~ n * a , 5(0) = 0; u = -1 + (1 - 2,tf/\ (4.2) 
(1 + it) 

which determines the scale for the solution. This scale together with the scaling 
invariance (|3.1|) . motivates the following expansion for the stretching boundary region 
in the vicinity of the end point x = 1 

u(x,t) = f(t)v(r,,t), V= J'* f{t) = l-{l-'St) 1 '\ (4.3) 

Note that f(t) = t+0(t 2 ) as t — > so an expansion of v(r], t) in powers of i corresponds 
at lowest order to an expansion in small f(t) and matches to the outer region exactly. 
Employing variables (|4.3[) together with the expansion 



k 

v( V ,t) = J2f n (t>n(v)+ £ (e 1/2 f(t) 1/4 )~v*(ri) (4.4) 

n=0 fc=l,fc^4p,p6N 



e 2 ' 



for the solution gives a sequence of problems to be solved for VkOq), k = 0, 1, 2, . . .. 

4 

The ©(e 1 / 2 /^) 1 / 4 ) component of (|4.4|) is the first correction to the 2r~ 1 u rrr term 
which appears, in the radial case N = 2, at a lower order due to the expansion not 
being centred on the origin. As can be seen from the equations below, when N = 1, 

all of the Vk with non-integer indexes may be chosen zero, so that the profiles Vk (77), 

4 4 

for k mod 4^0, play no role in the ID case. Equating powers of fit) 1 / 4 yields 

1, 77 > 0; (4.5a) 

(N - l)Gi (v (r)),vi(r)),vi(r]),V3(r))) + ~u 07) , V > 0; 

(4.5b) 

v 2n + 3v 2 =e 4 (N - l)G 2 fu (J7),wi(?7),ui(Y7),...,U7(77)J 

-3(v -r)^f+v^j+^v lv -2vi, t?>0; (4.5c) 

2(N-l)v 0rim , f?>0; (4.5d) 

(JV-l)Gi (yo(ji), , »7>0; (4.5e) 

(iV-l)G f (v (ri), vi (v),vl(v)) , V>0- (4-5f) 
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Fig. 4.1. Numerical solutions of equations 114.511 . On the left panel vq (solid 
curve), Vi (dashed curve) and vi (dotted curve) are displayed for N = 1 under clamped 
boundary conditions. On the right panel, vg (solid curve), vi (dashed curve) and vi 

4 2 

(dotted curve) are displayed for N = 2 under Navier boundary conditions. 



In the above, the functions G t represent lower order terms that only contribute when 
TV = 2. In what follows, we retain the first three non-zero terms of expansion (|4.4I) 

when TV — 1, and the first three terms (vn, vi and vi) when TV = 2. The above 

. 4 2 

equations are then solved together with boundary and far field behaviour 

(Clamped): Vj(0) = Vj v (0) = 0, Vjr,, Vj vvv -> 0, r; ->• oo, j = 0, 1, 2, -j, |. 

Uj(0) = Ujjjjj(O) = 0, v jn,Vjnnn -> 0, 7] -> OO, j = 0, 1, 2; 

(Navier) : 

W|(0) = V| w (0)-^ r) (0)=0, U|^,uj WJ7 ^-0, ?7->c5o, fc = L2 

_ (4-5g) 

The ODEs of (|4.5p are solved numerically as boundary value problems on an interval 
[0, L] with L taken to be sufficiently large so that their limiting behaviour for 77 — > 00 
is well manifested. Several profiles Vo(r)), Ui (77), 112(7?), V 1 (77), V2 (??), are displayed in 
Fig. 14.11 for both boundary conditions. 

In the ID strip case (N = 1), a solution valid for a; G (—1, 1) is obtained by super- 
imposing the left and right boundary phenomena and subtracting the extra far field 
solution to give a uniform approximation. This gives the small time approximate 
solution in ID case 

2 

u(x,t) = f(t)J2f n (t) 

71=0 

In Fig. 14.21 a comparison of the full numerical solution of (|1.5j) and the asymptotic 
solution (|4.6|) is displayed. Very good agreement is observed for small t. As t —> 1/3", 
fit) — > 0(1) indicating that the asymptotic solution (|4.6j) breaks down. Later in time, 
a new asymptotic regime based on small (i c — t) is entered. This touchdown regime 
is explored in iJS] 

In the unit disk case (TV = 2), the three leading terms in the asymptotic solution, 
each of which is displayed in Fig. 14.11 are 

u(r,t) = f(t)±(eVif(t)^) k v i y^)- (4-7) 



^1/2/(4)1/4 



1 - X 



£V2/(i)l/4 



-/(*) (4-6) 
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(a) e = 0.02 (b) e = 0.2 



Fig. 4.2. Comparison of full numerical solution ( solid line ) to l|1.5|l on ID strip 
with clamped boundary conitions to the asymptotic prediction ( dashed line ) of equation 
J4.6P . Left panel shows e < e c so that multiple touchdown points are present while panel 
(b) has e c < e < e* so that touchdown occurs at the origin. In both cases, solutions 
are increasing in time from top to bottom and good agreement between numerics and 
asymptotics is observed right up till the numerical solution enters the touchdown regime. 




(a) e = 0.02 (b) e = 0.2 



Fig. 4.3. Comparison of full radially symmetric numerical solution (solid line) 
to 11.51 on the unit disc with Navier boundary conditions to the asymptotic prediction 
(dashed line) of equation 1 14.711 . Left panel shows e < e c so that touchdown occurs on a 
ring of points while panel (b) has e c < e < £* so that touchdown occurs at the origin. 
In both cases, solutions are increasing in time from top to bottom and good agreement 
between numerics and asymptotics is observed right up till the numerical solution enters 
the touchdown regime. 

Fig. 14.31 displays a comparison of the full numerical solution to (|1.5[) and the asymp- 
totic solution (|4.7p . Good agreement is again observed for t small which breaks down 
as t — > 1/3 and the touchdown regime is entered. 

4.1.1. Estimation of Touchdown Points. To estimate the touchdown points 
of p.5[) , the critical points of the small t approximations (|4.6[) and (|4.7[) are examined. 
The first trough of the profile v(rj,t), defined in (|4.3p . serves as an estimator of the 
touchdown points and so its approximate value is determined asymptotically from 
(|4.4p . Minima of u are candidates for touchdown points with their location determined 
by the zeros of the derivative of u where, to leading order, the candidates for the 
touchdown points satisfy 
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e i/ay(t)i/4 J o ^1/2^)1/4 
1 - r 



£ V2/(i)l/4 



(N = 2) (4.8b) 



Note that when 1 + x c = O^ 1 / 2 /^) 1 / 4 ), 1 - x c = 0(1) and so for e 1 / 2 /^) 1 /* < 1, a 
zero of the left hand side of (|4.8aj) corresponds to the far field, i.e. flat, region of the 
right hand side of ()4.8a|) . In other words, for e 1 / 2 /(i) 1/ ' 4 <C 1 in the strip case N = 1, 
the two propagating regions do not interact directly and the critical points are the 
local maximums of the profile v(rj,t) inside each of the two regions. This assumption 
breaks down when e 1 / 2 /^) 1 / 4 = 0(1) as the two waves will superimpose to generate 
more complex solutions of (|4.8|) . 

The critical point inside each expanding region, rj c (t), satisfies 

(N = l) Vc(t)=Vo + f(t)vi+f 2 (t)V2 + --- /s , 

(TV = 2) Vc (t)= m + e 1 / 2 f(t) 1 / i r ]i +ef(t) 1 /\+--- 

where the corrections arc determined asymptotically from the condition v v (r] c (t), t) = 
0. In the N = 1 case, this provides the condition 

= v 0v (t] c ) + fv lr ,(r) c ) + f 2 v 2v (i]c) H 

= v 0ri (rio) + f[viVo vv (vo) + Vi v (rio)] 

v 2 

+ f 2 [v2 V (vo) + V2Vo, m (vo) + 77iwi w ( 7 to) + y«0wj(%)] H 

which gives the following definition for the corrections rjj, j = 0, 1, 2; 

vibrio) 



vonM = o, m 
-1 



m 



vo m (vo) 



v 2 

v2 V (vo) + m vim (m) + -^- v o,p P1 {vo) 



(4.9) 



A similar calculation can be performed for the N — 2 case and so the values of 
770 , 771 , 772 for N = 1 with clamped boundary conditions and 770,771,771 for N = 2 with 
Navier boundary conditions are found to be 

(Clamped): 770 = 3.7384, 771 = -0.6641, 772 = 0.1085 (4.10a) 

(Navier) : 770 = 2.8832, 771 = 0.3533, 771 = 0.9457. (4.10b) 

This now allows for the two critical points xf(t) in the strip case N = 1 and the 
ring of touchdown points r c (t) in the radially symmetric unit disc case N = 2 to be 
specified as 

IV = 1 xt(t) = ±[l-e 1 / 2 f(t) 1 ^[ m + f(t) Vl +f 2 (t) m ]\, (4.11a) 
N = 2 r c = 1 - e 1/2 /(ic) 1/4 /7o - e/(tc) V V - £ 3 / 2 /(*c) 3/ V + • • ■ (4.11b) 
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Fig. 4.4. Touchdown location for 11.51 1 from full numerics ( solid line ), compared 
with asymptotic formula 114.111 1 with t c from full numerics ( dashed line ) and asymptotic 
formula lj4.11| l with t c = 1/3 ( dotted line ). Left figure case N = 1 with clamped 
boundary conditions, right figure case N = 2 with Navier boundary conditions. 



Note that the approximation for the touchdown locations requires t c , the touchdown 
time of (|1.5[) . As observed in Fig. 14.41 asymptotic formula (|4.11[) captures the location 
of touchdown very well, particularly when s -C e c - As e — > e~ , the approximation 
breaks as the left and right boundary effects are superimposing and so the touchdown 
points are no longer simply the minima of the isolated profile v(y, t).ali 

5. Touchdown regime. To establish a blow up profile in the touchdown regime, 
the techniques of |22j are employed. The correct similarity variables are investigated 
by initially rescaling equation (jl.5ap with 

u = — 1 + Uu(x, t), t = Tt, x = Lx 

which results in 

U . U 



A balance of all terms suggests scaling with L ~ y 1 / 4 and U ~ T 1 / 3 and so an 
appropriate self similar solution would be of form 

1/3, 



■i+R{ty> A v 



where R(t) is the quenching rate of the solution. In general, rigorous determination 
of R(t) is a difficult problem and so we make reasonable guesses and investigate their 
validity with numerical calculations. This approach is not definitive, however, as the 
case of blow-up in critical NLS[5], whereby the rate has been found to satisfy the so 
called loglog law 

2n(t c - t) 

indicates. Numerical verification of this rate law would require accurate solutions for 
almost surely unobtainably small values of (t c — t). As such, the evidence presented 
here for self-similar quenching awaits rigorous verification. The case of quenching 
solutions in the strip and unit disc geometries are treated separately and both appear 
to be self-similar in nature. 
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5.1. Touchdown solutions in ID. In similarity variables 

u(x,t) = -l+(t c ~t) 1 / 3 v( V , S ), n= e y 2 X ~*' /4 , « = -l0g(t c -t) (5.1) 

(ll.5a|) is transformed to 



77 v 1 
— Vn H — 
4 ' 3 



v v+o-7a> (?7,s) £ Kx M+. (5.2) 



Far field and initial conditions for v(rj, s) are now discussed. The behaviour of 17(77, s) 
for 77 — > ±00 corresponds to a solution of u(x,t) for x 7^ x c as t — > t~ . Assuming a 
localized quenching solution at x = x c , it can be expected that u t = 0(1) in a region 
away from x c as t — > t~ . Now, 



« t = (tc-t)- 2/3 

and so the condition that ut = 0(1) implies that 



77 v 

— Vn 

4 11 3 



(5.3) 



Vs + \v n - V - = 0{{t c ~t) 2 /*), t^t- (5.4) 

For a fixed x ^ x c , the limit t — > t J corresponds to |?y| — > 00 and so (|5.4[) augments 
(|53]) to establish 

v s = -v mm - -v v + - - -j , (77, s) 6 K x M+; (5.5a) 
u 77 

u s = - - -v ri , 77 ±00; (5.5b) 

A key step is to determine the limiting behavior of solutions to (|5.5[) for any fixed 77 
as s — > 00. One obvious candidate for an equilibrium state is the constant v = 3 1 / 3 . 
An analysis of its stability leads one to consider the eigenvalue problem 



for m — 2. The spectrum of the operator C m in the weighted space L^(M.) where 
p = e~ a W , with a positive ( c.f. [HE]), is 

a(£ m ) = ^ k = 1 - ^; fc = 0,l,2,...J (5.7) 

and so there are two linearly unstable modes associated with this equilibrium, po = 1 
and fii = 1 — 1/m for m > 2. The instability associated with the mode fj,o = 1 
is generated by the invariance of the touchdown time t c and is therefore not a true 
instability. The instability associated with the fii = 1 — 1/m mode represents a 
true instability when m > 2. The presence of this positive eigenvalue indicates that 
equation (|5.5p does not satisfy v(r],s) — > 3 1 / 3 for fixed 77 as s — > 00, thus we seek 
equilibrium solutions of the following nonlinear problem 



77 _ v 1 

v vvvv + 7 W ») ~ 3 + =2 = °) ~°° < ^ < 00; (5.8a) 



I - = 0, 77 -> ±00. (5.8b) 
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and investigate the multiplicity and stability of its solution. 

The robin condition of j5J3b| suggests that (|5.8p admits a far field series solution of 



forr 



Cn\V 



4/3-4n 



\rj\ — > ±00. 



(5.9) 



n=0 



Here the constants c„ = c„(co) are functions of the parameter cq for n > 1 and can be 
determined by lengthy but straightforward manipulations, e.g. c\ = 40co/81 + c^ 2 . 
The parameter cq plays the role of a nonlinear eigenvalue and it is expected that 
(|5.8|) will have solutions for isolated values only. Taking the limit t — > t~ for fixed 
x 7^ x c corresponds to the limit |?7| — ¥ 00 and therefore, in physical co-ordinates, the 
touchdown profile is expected to satisfy 



w(x, t) 



T + c 



rl/2 



4/8 



+ Cx{t c -t) 



rl/2 



-8/3 



+ ■ 



a.s 



t -> t~ (5.10) 



Additional boundary conditions are now obtained for (|5.8j) by suppressing exponen- 
tially growing modes of the linearization of ()5.8|) about v for large 77. To analyze 
linearized perturbations of (|5.8[) about v p , set v — v p + aw where a < 1 to arrive at 
the equation 



1 V 

^rfrfrfr/ T ^ 



^-2^ = 0, 



For large |??|, a WKB anzatz solution of form 

" 1 00 



w 



exp 



fc=0 



-co < 77 < OO. 



c 

n = 

v 



(5.11) 



for v <C 1 and <5 



.4/3 



produces the leading order equation 
9oc + l9o( = 



which admits three exponential solutions 

30,(0 = -3|C| 4/3 2- 8/3 exp 



2nij 



J =0,1,2. 



The terms exp(<?o.j) for j = 1,2 are growing as rj — > ±00 and need to be suppressed in 
the solution of (|5.8j) . The mode corresponding to g' = is w = ?y 4 / 3 which represents 
an arbitrary change in the value of cq. At the following order exp(<7i(£)) = ?y -10 / 9 
which now gives the following full specification for v(r]) 



0, 



-00 < 77 < 00. 



"' 4 " 3 i 

DO 

71 ~ c„77 4 / 3 - 4 " + ClTyr 10 / 9 exp [-3M 4 / 3 2- 8 / 3 



n=0 



(5.12a) 
77 -> ±00 (5.12b) 



Extracting information from (|5.12j) is analytically challenging as it involves solving 
a fourth order, nonlinear, non constant coefficient and non-variational differential 
equation. This motivates the use of numerical techniques to analyze the multiplicity 
and stability of solutions to (|5.12[) . 
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(a) Profile vi(n) (b) Profile V2(rj) 



Fig. 5.1. Plots of two self-similar profiles v\(r]) and V2{rf) satisfying {5TT2J- The 
dotted curves represent the far field behaviour Vj (n) ~ Cq |t;| 4 / 3 as \r/\ — > oo. The values 

Cq 1 ' = 0.906 and c^ = 0.1047 were determined numerically. Note that V2 has a small 
dimple at the origin indicating three critical points. 



5.1.1. Numerical and stability analysis. This section deals with the numer- 
ical determination and linear stability of solutions to (|5.12[) . Related similarity ODES 
have been solved by several authors in the context of pinch off dynamics for thin films 
[2"5I [53] and a framework for their solution is well established. Equation (|5.12al) is 
solved by first applying a centered difference discretization scheme to the derivative 
terms on a uniform grid of [—L,L]. The Robin condition (|5.8b[) is discrctized and 
applied to remove the ghost points from both end points and thus effectively yields 
four boundary conditions for the system. The application of the Robin condition at 
two nodal points enforces the far field behaviour v ~ co|?7| 4 / 3 and also eliminates 
exponentially growing terms. 

This discretization leads to a large system of nonlinear equations to be solved via 
a relaxed Newton's Method [JJ. The iterations are initialized with a solution of the 
reduced equation 

\»n ~ I + ^2 = 5 = + c o > (5.13) 

over a wide range of positive Co until convergence is achieved. This initial guess has 
the advantage of satisfying the far field behaviour exactly for a given Co and also being 
smooth at the origin. The size L of the system is taken to be sufficiently large so that 
the far field behaviour is well manifested. After seeking convergence over a wide range 
of parameters cq, exactly two solutions to (|5.12j) . denoted vi{rj), v 2 {r]) were found as 
shown in Fig. 15.11 This solution multiplicity appears to be qualitatively similar in 
character to that observed [2] [IT] [12] in the self-similar blow up of fourth order PDEs 
with power law nonlinearity. To address the question of the existence of a stable 
self-similar quenching profile for (|1.5j) . the linear stability of vi{rj) and V2(t]) is now 
analyzed by setting v = v(rf) + 4>(r])e^ s for <j> <C 1 in (|5.5[) to arrive at the eigenvalue 
problem 

fxcf) = -<t> mm - 7 -t<f> v + Q + 0, -oo < V < oo; (5.14a) 
/U<?!> = - - r/^±oo. (5.14b) 
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Mo 


Mi 


M2 


M3 


M4 


M5 


Me 


M7 


Vl 


1.0003 


0.2499 


-0.1369 


-0.4328 


-0.6089 


-0.8431 


-1.1431 


-1.4251 


V2 


1.0000 


0.7740 


0.5347 


0.2499 


-0.0828 


-0.4464 


-0.8269 


-1.2151 



Table 5.1 

The first eight numerically obtained eigenvalues of l|5.14[l for the two profiles vi (n) 
and V2(rj). The value of L = 50 and a uniform discretization with N = 1000 grid points 
were used. 




Fig. 5.2. Comparison of full numerical solutions (dashed lines) of 11.51 to the 
stable self similar profile (solid line) vi(n) for various t approaching touchdown at t c . 
This is for N = 1 and clamped boundary conditions. 



Apart from the following two modes associated with translation in touchdown time t c 
and location x c 

V Tj _ 1 

Mo = 1, ( / , o=g-^-w, ) ; Mi = ^7 4>i =v„, (5.15) 

the spectra of (|5.14j) must in general be determined numerically by reducing, via 
discretization, (|5 . 14[) to a linear system C^4> = and then seeking /i such that det = 
0. The eigenvalues appear to be purely real and the largest eight numerically obtained 
eigenvalues associated with each of the two profiles v\ (rj) and v?, (77) are displayed in 
Table. 15.11 In the spectra associated with each profile, the two eigenvalues identified 
in (|5.15p are present. Ignoring these particular values, it is observed that the spectrum 
associated to v\ is strictly negative, while the spectrum associated with contains 
two positive eigenvalues. 

This suggests that the profile v\(r\) is a stable self-similar quenching profile for (|1.5[) 
and indeed, in Fig. l5.2[ convergence of the full numerical solution to ^1(77) is observed 
as t — > t~ for the case of touchdown at and away from the origin. 

5.2. Radially symmetric quenching solutions in 2D. Self similar quenching 
profiles of the MEMS problem (|1.5[) are now considered in two spatial dimensions. For 
radially symmetric solutions on the unit disc, the cases of touchdown at and away 
from the origin are treated separately. The variables 

U(X, *) 1 + (*e - ^MV), V = i/ 2r /_,U/4 ' ( 5J6 ) 
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(a) Profile v\ (n) 



(b) Profile V2(r]) 



Fig. 5.3. Plots of two self- similar profiles vi(r)) and t>2(p) satisfying I I5.19H . The 
dotted curves represent the asymptotic far field behaviour Vj(p) ~ Cq p 4 / 3 as p — > oo. 



The values c, 



(l) 




0.7265 and c, 



(2) 




0.0966 were determined numerically. Note that V2 



has a small dimple at the origin indicating two critical points including that at p = 0. 



which assume touchdown at the origin, transform (|1.5[) to 



1 



1 



7] e 



(5.17) 



which is a partial differential equation for the self similar quenching profile. The 
question of existence, multiplicity and stability of solutions to (|5.17|) appears to be an 
open question in spatial dimensions N > 2 (c.f. |12j). If a radially symmetric solution 
v{rf) = v(\r)\) is presumed, then (|5 . 1 7[) reduces to 



P 



~ A pv 



V 1 

3 v z 



(5.18) 



where p = \r]\. A far field analysis similar to that which led to (|5.12[) can be applied 
to (|5.17[) to establish boundary conditions which imply algebraic growth with expo- 
nentially growing terms suppressed at infinity. After algebra the full problem for the 
radially symmetric self-similar quenching profile in dimension N = 2 is 



//// . 

v + 



v(p) 



coP 



P 

4/3 



// 1 

P 2 



1 



1 



p 3 4 

Cp -16/9 



pv 



cxp 



-3p 4 / 3 2- 



= 



•8/3 



with symmetric conditions at the origin enforced by 

w'(0) = v"'(0) = 0. 



p > (5.19a) 
p^oo (5.19b) 

(5.19c) 



This nonlinear equation is solved numerically by first discretizing ()5. 19a[) on [0, L] 
for L large, applying far field behaviour (|5.19b[) as a Robin condition ( c.f. (|5.8b[) ) 
at consecutive endpoints followed by Newton iterations initialized with (|5 . 13|) . The 
iterations are initialized over a wide range of the parameter cq, with convergence 
observed for two isolated values = 0.7265 and Cq 2 ' = 0.0966. The two associated 



profiles are displayed in Fig. 15.31 

As in the N = 1 case, the second profile V2{ff) has a dimple at the origin and, as 
illustrated in Fig. 5.4(a)| full numerical solutions of (|1.5|) are observed to converge to 
the monotonic self-similar profile vi{rf). 
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Fig. 5.4. Convergence of radially symmetric solutions (dashed & dotted) of 11.51 
to self-similar profiles. Left: Touchdown is at the origin and convergence is observed 
to the monotone profile vi(n) (solid) solving 15.191 1. Right: For this case, touchdown 
is away from the origin and so convergence is to the monotone profile vi(t}) solving 
l|5.12( . as in the ID strip case. In both figures, the dotted curve represents the solution 
for smallest (t c — i). 



For touchdown away from the origin in the radially symmetric unit disc case, the 
self-similar quenching profile appears to be the same as that obtained for the N = 1 
case. Indeed, the appropriate similarity variables are 

u(x, *) 1 + (*e - t) 1/3 v(v), V = £ l/2( fe ~_^)l/4 - ( 5 - 20 ) 

These variables rescale the biharmonic term as follows 

—e 2 A 2 u -+ -(t c t)- 2/3 + 0{e^l\t c - 

and so in the limit as t — > t c , the u ww term is dominant. This results in a self-similar 
profile which satisfies 



Tj _ v 1 

t -v„ - - + -rj = 0, -oo < i] < oo. (5.21a) 



v 



CnV i/3 ~ in + C\ V \- W / 9 cxp [-3M 4 / 3 2- 8 / 3 ] , 77 -> ±oo (5.21b) 



n=0 



as derived for the ID case in (|5.12[) . Consequently, quenching solutions away from the 
origin in the radially symmetric unit disc case are expected to converge to the self- 
similar quenching profile of the N = 1 case, as confirmed by the numerical simulations 



displayed in Fig. 5.4(b) 



6. Discussion. Quenching solutions of a fourth order parabolic differential equa- 
tion with a singular nonlinearity have been analyzed for a ID strip and under radial 
symmetry on the unit disc with both clamped and Navier boundary conditions. In 
contrast to its second order equivalent, the fourth order PDE can quench at multiple 
points away from the origin. More precisely, in the case N = 1, we have shown that 
the PDE can quench at two distinct points symmetric about the origin, while in the 
radially symmetric unit disc case, it can quench on an inner circle of finite radius. In 
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each case, the location of the quenching set was predicted by means of an asymptotic 
expansion whose accuracy is verified using very accurate adaptive numerical methods. 

In the limit as t — > t c , where t c is the quenching time, the behaviour of (|1.5[) was 
shown to be self-similar in nature. This is again in contrast with the second order 
equivalent of ()1.5|) . The self-similar profile itself was obtained numerically and its 
limiting behaviour for £ — ^ i c is given by 



where c = 0.7265 for N = 2 and touchdown at the origin and Co = 0.9060 in the 
other cases. 

There are many interesting questions which stem from this study. In the case of the 
unit disc, it may be possible for the dynamics of (jl.5l) to break the radial symmetry of 
the quenching set. All the simulations presented here were initialized with u(x, 0) = 0. 
Adding random noise to the initial condition breaks the left-right symmetry for N = 1 
and rotational symmetry for N = 2. The symmetry breaking can be amplified by the 
dynamics of the PDE. 

In such a scenario, the ring would most likely be split up into a collection of points 
whose arrangement would need to be determined. The prediction of the quenching 
set of (|1.5p for larger classes of 2D geometries is another interesting open problem. 
For regular geometries, it may be that the number of axes of symmetry determines 
the quenching set but for irregular domains, it is not clear that the touchdown lo- 
cations can be determined by simple geometric considerations. This question may 
be amenable to perturbation analysis, for example an almost circular domain whose 
boundary is r = 1 + Sf(6) for some i < 1 and f{ff) a 27T periodic function. 

A robust method for solving (|1.5jl of a large class of 2D geometries would be an essen- 
tial complement to any analytical investigation of the above questions. In particular, 
a meshlcss method might be well suited to handle the highly non-uniform grids needed 
to resolve the dynamics of (|1.5|) very close to touchdown [8]. 

The treatment of these open issues is beyond the scope of this manuscript and will 
be left for future investigation. 

7. Acknowledgements. A.E.L is very grateful to M.J. Ward for many useful 
discussions. 

Appendix A. Spatial Discretization. 

For discretization in space, a collocation method based on piecewise 7th-order polyno- 
mial interpolation is employed. On the interval x £ [Xi(t), Xi + i(t)] for i = 0, 1, . . . , N, 
the solution u(x,t) is written as 




3 



u(x, t) 



Y^[nt\t)L M+u^ 1 {t)L 1 , k {s i )]H^ 



(1.1a) 



k=0 



where 



x-Xj(t) 
Hi(t) 



e [0,1], 



Hi(t) = X i+1 {t)-Xi{t), 




dx~* uM 



J x=Xi(i) 



(1.1b) 
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and the Lq j(s), L 1 j(s) for j = 0, 1, 2, 3 are the shape functions 

L 0j0 ( s ) = (20s 3 + 10s 2 + 4s + l)(s - l) 4 , Lo.i(s) = s(10s 2 + 4s + l)(s - I) 4 

„2 ,,3 



LoA*) = y(4s + l)(s-l) 4 , 
Li, (s) = — s 4 (20s 3 — 70s 2 + i 



35), 



io, 3 ( S ) = y( S -l) 4 , 
^(s) = s 4 (s - l)(10s 2 - 24s + 15), 



Li, 2 (s) = - y (s-l) 2 (4s-5) 



They satisfy 



d p 

— — Li k(si) 
dxP ' v ' 



x=Xj (t) 



L li3 (s) = -(s-l) 3 . 



1 if i = j and k = p 
otherwise 



(1.1c) 
(Lid) 



so that the unknown coefficients iS (t) are the values of u and its first three spatial 
derivatives at the nodal points x — By construction, these are continuous at 

the nodal points. 

The dynamics of the u\ (t) is obtained by substituting expansion (|l.la[) into the 
PDE, using the following expressions for the temporal and spatial derivatives of u. 

3 r 



d_ 

di 



i(x,t) = j2 

fe=0 
3 

i(x,t)=j2 

k=0 L 

dH, 



u ( k \t)^L As l )+u^ 1 (t)j-L 1 , k (s l ) 



d-> 



H 



k-j 



(1.2a) 



± u f(t)L 0tk (s i ) + jU < S 1 {t)L lik (s i ) 



H 



fc=i 

. . dXi dH,, 
u x (x,t) —j— + Sf 



kH 



fc-i 



(1.2b) 



dt 



Navier and clamped boundary conditions can be applied at both endpoints by choosing 



l N+l 



(2) _ ,.(0) 

U N+1 



4 AH 
,(2) 



U N+1 = 



U N+1 = 



,(!) 



(Clamped) = 
(Navier) Uq = u Q 
in the strip f2 = [—1,1] case and 

(Clamped) Uq 
(Navier) Uq = u 1 ^' = u^' +1 = 
in the unit disc f2 = {x 2 + y 2 < 1} case. 

and the remaining equations are obtained by writing the discretized PDE at the Gauss 
points 



(3) 
U Q 


(0) 

— U N + 1 ~ 


u (1) 

Ll N+l 


(0) 
N+l 


-u i2) + 


U N+1 



= 
= 



Pl = 2 



1 V525 + 7CK/30 



70 ' 92 2 70 

on each interval [Xi(t), Xi + i(t)], for i = 0, . . . , N. This provides 4(iV + 1) equations, 
which together with the four boundary conditions, are integrated in time to obtain 



1 ^525 - 70\/30 



P3 = 1 - P2, pi = 1 - Pi 



the 4(7V + 2) unknown nodal values u\ k \t) for k = 0, 1,2,3 and 



0. 
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